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Abstract 

The asymptotic travelling wave solution of the KdV-Burgers equation driven by 
the long scale periodic driver is constructed. The solution represents a shock-train 
in which the quasi-periodic sequence of dispersive shocks or soliton chains is inter- 
spersed by smoothly varying regions. It is shown that the periodic solution which 
has the spatial driver period undergoes period doublings as the governing parame- 
ter changes. Two types of chaotic behavior are considered. The first type is a weak 
chaos, where only a small chaotic deviation from the periodic solution occurs. The 
second type corresponds to the developed chaos where the solution "ignores" the 
driver period and represents a random sequence of uncorrelated shocks. In the case 
of weak chaos the shock coordinate being repeatedly mapped over the driver period 
moves on a chaotic attractor, while in the case of developed chaos it moves on a 
repellor. 

Both solutions depend on a parameter indicating the reference shock position in 
the shock-train. The structure of a one dimensional set to which this parameter 
belongs is investigated. This set contains measure one intervals around the fixed 
points in the case of periodic or weakly chaotic solutions and it becomes a fractal 
in the case of strong chaos. The capacity dimension of this set is calculated. 
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1 Introduction 



Both the Korteweg de Vries (KdV) and the Burgers equations are widely used 
in the study of weakly nonlinear wave phenomena in a variety of continu- 
ous media. These equations have the same type of nonlinearity and represent 
dispersive and viscous system, respectively. In many realistic cases, however, 
both viscosity and dispersion are equally important (see eg. [1-3] and references 
therein) and a generalization of these two equations known as the Korteweg- 
de Vries -Burgers equation (KdVB) is very useful approximation to study the 
wave dynamics. 

For wave motion to appear the equilibrium state of the system must be dis- 
turbed. Traditionally, this is attributed to initial conditions for the undriven 
equations. As far as the KdV and the Burgers equations are concerned, the 
time asymptotic behavior of the solution can often be described in terms of the 
dynamical evolution of fundamental solutions of these two equations, which 
are the cnoidal waves (the solitons in the limiting case) and the shocks for the 
KdV and for the Burgers equation, respectively (see e.g. [4,5]). 

Different possibilities arise in the case of a continuously acting driver, which 
can be of an external or a growth rate type. The solutions of the undriven KdV 
and Burgers equations can also play a certain role when the driver is switched 
on. The driven equation itself, however, can possess a fundamental solution 
which is completely different from any solution of the undriven equation. Such 
a solution, called shock-train, has been proposed [6,7] in the framework of the 
so called derivative nonlinear Schroedinger-Burgers (DNLSB) equation, which 
is the viscous version of the DNLS-equation, introduced by Rogister, [8]. The 
idea behind this solution can be described as follows. In the case when both the 
dispersion and the viscosity are small, the solution is defined throughout most 
of the driver period by the balance between the nonlinear term and the driver. 
Since the nonlinearity of the equation provides more than one equilibrium 
state (two in case of the KdVB and three in case of the DNLSB equation, 
depending on the order of the nonlinear term), the overall solution can in 
principle be composed of the pieces of these different equilibria, connected 
by thin regions (shocks) where high derivative terms take part in the total 
balance. In terms of the ordinary differential equation (ODE) to which the 
original evolution equation is reduced when one is looking for the travelling 
wave solution, these shocks correspond to the transition between the unstable 
and stable types of equilibrium. Repeating this in the subsequent periods of 
the driver, we obtain a sequence of shocks, the shock train. Depending on the 
values of dispersion and viscosity each shock can be very oscillatory and can 
look more like a chain of solitons than a shock. The major difficulty here is 
to define the shock locations relative to the driver phase. This problem can 
be put in the following context. Suppose we know the location of the shock, 
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if any, within a given driver period which is equivalent to the specification 
of the boundary condition for the ODE. The question then is, where is the 
shock located within the neighbouring period? Provided that we can answer 
this question, the solution can in principle be calculated everywhere. 

For the Burgers equation the answer is very simple: each shock must be at 
the same position relative to the driver phase, which yields a strictly periodic 
solution having the driver period [9]. This can be understood from Poincare- 
Bendixson theorem in terms of the insufficient dimensionality of the related 
ODE. Indeed, the dynamical system that describes this travelling wave so- 
lution is two dimensional, and according to Poincare-Bendixson theorem its 
attractors are limited by the fixed points and the limit circles. More extensive 
discussion of this issue and the relation to another type of solution that consists 
of choosing a stable branch of driver-nonlinearity equilibrium, can be found in 
[10]. The latter solution, containing no shocks at all has been constructed in 
[11] in the framework of the periodically driven KdV equation. In the present 
paper we demonstrate the possibility of a completely different response to an 
external force, giving rise to the formation of a shock-train. Moreover, this re- 
sponse seems to be more persistent, due to its ability to dissipate the pumped 
energy even at low viscosity, since shocks are involved. 

Like other methods of description of complicated nonlinear systems, the shock- 
train technique outlined above can in principle produce physically irrelevant 
solutions. In particular, the shock train method reduces the continuous so- 
lution to a Poincare map of the shock coordinate. This map is no longer 
constrained by the Poincare-Bendixson theorem (even when the original dy- 
namical system is) and in principle can falsify spatially chaotic behavior. An- 
other danger is the possibility of the artificial temporal chaos caused by an 
insufficient number of modes in the computation of the shock train dynam- 
ics. For example such chaos has been observed during the numerical study 
of the unstably driven Burgers equation [10]. This observation is interesting 
in view of the remarkable and, after Lorenz's [13] famous paper widely ex- 
ploited fact, that finite mode representation of continuous media can reveal 
chaotic behavior even for a very modest number of modes. In the case of 
shock-train dynamics the number of modes must be rather large to properly 
describe a shock structure at viscous and dispersive scales. The truncation 
virtually introduces the strong damping at the short scales which facilitates 
the appearance of stochastic attractor and rules out the shock-train, whereas 
the full solution (i.e. calculated with sufficiently large number of modes) shows 
that only regular shock-train dynamics occurs. Coming back to the possibility 
of spatial chaos, we note that the KdVB travelling waves are governed by the 
driven-damped nonlinear oscillator equation (see eq.(4) in the next section) 
which has a three dimensional phase space, and hence is not constrained by 
the Poincare-Benidixson theorem in contrast to the Burgers equation. 
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A similar problem of spatial chaotization of periodically driven Alfven waves 
has been studied numerically in the framework of the DNLSB by Hada et al. 
[14]. In the present paper we turn to the KdVB equation which is a subset of 
the DNLSB and in its context describes weakly nonlinear oblique MHD waves. 
We also restrict ourselves to the case of small viscosity and dispersion, which 
is interesting for two reasons. First, this case corresponds to a quite typical 
situation in study of turbulence when the energy is pumped into a system at 
long scales, then transferred over an inertial range to short scales via non- 
linear interaction of modes, and finally is dissipated there due to viscosity or 
any other short scale dissipation mechanism. Second, this case is difficult for 
numerical integration since different spatial scales appear. A direct manifes- 
tation of spatial chaos in the time dependent numerical study would formally 
require an infinitely large computation box or infinite number of modes. 

The plan of the paper is as follows. In section 2 the shock-train solution 
to the KdVB equation within a driver period is obtained. In section 3 the 
solution from two neighbouring periods are matched. In section 4 the mapping 
of the shock coordinate at the driver period is introduced. In section 5 a route 
to chaos through the period doubling bifurcations is obtained as well as a 
developed chaos in the mapping of the shock coordinate is considered. Section 
6 deals with numerical integration of the driven KdVB equation and compares 
the results. 



2 Solution within one driver period 

2.1 Overall behaviour 

The KdVB equation that describes nonlinear waves, driven by a source moving 
with the constant velocity can be written in the following form 



We have transformed the coordinate system to the reference frame where a 
periodic driver 




(1) 



Q' = dQ/dx 



is at rest; u and x are normalized in such a way, that 



Q(x + 2ir) =Q(x), 
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and the function Q is of order unity whereas the dispersion coefficient A 2 and 
the viscosity fi are both assumed to be small. We removed the linear term 
from the l.h.s. of eq.(l) by shifting u — > u — c/2. From the physical point 
of view this choice of x and u variables corresponds to the case when the 
driver moves at the speed of a linear nondispersive undamped (A, /i — > 0) 
wave and, therefore, is in resonance with it. However, this does not mean 
the loss of generality, unless u is constrained by additional conditions like 
x — > ±00 asymptotics. Note, that this particular constraint is relevant rather 
for the unforced case or for the case of localized forcing. In the case of periodic 
forcing, considered here a more suitable constraint is the mass invariant 
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that can be introduced for the 2nn- periodic (n is integer) solution and is 
affected by the shift in u. In the above context Eq.(l) remains generic for all 
steady drivers unless u in Eq.(2)is specified. Below we specify explicitely a 
class of particular solutions of Eq.(l) considered throughout the most of this 
paper and we shall briefly return to their dependence upon the mass invariant 
(2) in sec. 6, where these solutions are verified by the numerical integration of 
Eq.(l). 

The nondispersive shock-train solutions do not allow the spatial chaotization, 
as it was stressed in the preceding section. We therefore consider the case of 
relatively strong dispersion 

At « A < 1 (3) 



when the developed oscillatory structures must appear nearby the shocks. At 
the same time, intending to obtain a solution with distinguished shocks in the 
shock-train we assume that those oscillations are appreciably damped between 
two shocks. The spatial damping rate of these oscillations is defined by fi/X 2 
and we assume that 

P = /i/A 2 > 1. (4) 



Note, that under the conditions (3,4) an efficient method for Eq.(l) would 
be that one, developed by Whitham [15], and its generalization to the viscous 
driven case seems to be tractable. However, seeking the travelling wave solution 
of Eq.(l), which moves at the driver speed, i.e. assuming 

du/dt = 0, 
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we immediately arrive at the ODE which is the equation of the driven-damped 
nonlinear oscillator 

u 2 - \ 2 u xx - fj,u x = Q(x). (5) 



Our analytical treatment in what follows will concern only (5), but we shall 
return to Eq.(l) in sec. 6. 

Interesting in the shock-train solutions of Eqs.(l) and (5), we choose the in- 
tegration constant in the r.h.s. of (5), as 

Q{x) = Q {x) + Xd + fiC 2 + ■ ■ ■ , (6) 



where Cj are constants of order unity. According to the general shock-train 
method [9,10] these constants should be chosen in such a way, that Qo(x) > 
and Qo, Q' go to zero simultaneously. This admits a pair of inviscid solutions 
±VQo (to the leading order in A and fj) and the smooth transition between 
them at points x (mod 2n) where 

Qo(xo) = Q'o(xo) = 0, 



whereas each shock in the shocktrain is supposed to be somewhere between 
two such points. 

To begin, we introduce the fast variable 

Z = \9i(x)+92 + ---, (7) 

where gi(x) are functions to be found from the further analysis. We seek the 
solution of (5) in the following form 

u — u (£, x) + Aiii -I . (8) 



Then, the first two orders of approximation yield 



2 _ , 2 d 2 Mo 
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Qo 



(9) 



2u u 1 -g 1 -^ r = 2g 1 g 2 — + 



duo 



^ d_ , ( du 
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(10) 
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The solution of Eq.(9) is a cnoidal wave (e.g. [4]), which depends on i as on a 
parameter. For convenience of the further analysis, we first represent it in the 
following form 



«o(£,a0 = 2JQ 



c-i + (Co - c-i)sn 2 [ dy (Cl 3 C " l} ^fc 



where ( n are the roots of the cubic equation 

C 3 -tC + t = o, 



Cn = sin 



1 2tto 
- arcsm a H — — 

o o 



(12) 
(13) 



where n — 0, ±1, and fc denotes the modulus of the elliptic function in Eq.(ll) 

Co - C-i 



k 



Ci-C- 



(14) 



We have introduced the function | a(x) |< 1 or, equivalently, < k(x) < 1 
which is related to a(x) by Eqs. (12-14), to be determined as a constant of 
integration of Eq.(9). In order to avoid the appearance of secular terms in the 
r.h.s. of Eq.(10), we specify gi(x) in such a way that the period of the solution 
Uq(^,x) in £ is independent of x. Then, the x-dependence of duo/dC, in the 
r.h.s. of Eq.(10) will not produce a nonperiodic in £ term. Since 



G-C- 



2 



(15) 



we must set in Eq.(ll) 



9i = 



2iK(P + A;' 4 )3 



(16) 



where k' 2 = 1 — k 2 and K is the complete elliptic integral of the first kind. 
Eq.(ll) then transforms to the following 



1 Qo(x) 
k 2 + k' A 



-1-k 2 + 3rsn 2 (2K(fc)f,fc)) 



(17) 



To determine the function k(x) we multiply Eq.(10) by du /d£, and integrate 
over one period regarding x as a parameter. This yields the following simple 
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equation for J(x) 
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where 



=<*/(§*)"*■ (19) 



si 



Substituting Eqs.(16) and (17) J(x) rewrites 



J{x) = 



6V2Q 5 /4 (x) 
5(k 2 + fc /4 )V4 



2E(fc) 



k >2 (l + k' 2 ^ 
k 2 + k' 4 



-K(fc) 



(20) 



where E is the complete elliptic integral of the second kind. Eqs.(18) and 
(20) implicitly define the function k(x) in the solution (17). Note, that the 
integration constants in the solutions of Eqs.(16) and (18) are not specified 
yet. It is convenient to attribute one of them to the shock coordinate x* so 
that from Eqs.(7) and (16) we obtain 

c^9i(x)_ 1 J Qldx 

? A 23/ 2 A7 (k 2 + k'^K(k)' { } 



The function k(x) is defined by the following equation, which results from 
Eqs.(18) and (20) 

Ql(x)QoH<) 
2(k' 2 + k^ 



,-/3(x-x'») 



(22) 



where x'^ represents an arbitrary constant in the solution of Eq.(18), which we 
have specified in such a way that k(x'^) = 1, whereas < 1 (see Fig.l), 

and is the nearest to the point x'^ root of the function duo/d^. The solution 
(17) with k(x) defined by Eq.(22) cannot be continued to the left from x — x'^, 
since k exceeds unity in this region and (17) should be replaced there by a 
different solution (see Eq.(38) below). However, as we will see Eq.(17) yields 
the correct asymptotics in the region x[ — x ^ A, provided that k is set to 
unity at x < x'^. Therefore, we can define k(x) for x G (x — 2tt, x ] as follows 



k = I 



k(x), x'x < x < x 

1, Xq — 2-k < x < x'*, 



(23) 
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Fig. 1. An arrangement of the points x* and x'* in the zero order solution (17). Three 
solitons marked by digits fill up an interval between the shock front and the 'turning 

1/2 

point' xq. Inviscid nondispersive solutions ±Q are plotted with thin lines. 

where k(x) is the solution to Eq.(22) at x > x[. Indeed, since sn in Eq.(17) can 
be replaced in the region x < re* by tanh, (by virtue of A;' 2 (re*) ~ ^/i/A <C 1, 
see Eq.(36) below) the solution uq tends to the unstable (for growing x) branch 
of the inviscid nondispersive solution of Eq. (5) u ~ \[Q- At x > x* it oscillates 
around the stable branch u ~ —\fQ, Fig.l. Not too far from the points x* ~ x'^, 
i.e. at x—x'^ < 1/ [3, where k(x) is still close to unity, these oscillations virtually 
form a soliton chain attached to the shock front to the right from x' m . The 
number of these solitons is defined by the "optical length" between rr* and x . 
Beyond x$ the oscillatory behavior of the solution is changed for the monotonic 
one (see the next section). 



2.2 Resolving unsmooth region 



Although the solution (17,23) is smooth in the fast variable £, it is not smooth 

in the slow variable x at OC Ob ^ Ob ^ • The adiabatic approximation which 

allows the introduction of the fast and the slow variable breaks down there 
(the period of the solution in £ tends to infinity near the separatrix k — 1), 
although the unsmootheness of Uo(x) does not occur in the main order of 
approximation in A 1. A condition for the smooth matching of the solution 
at x = x'*, where the derivative of k(x) has a discontinuity, links the constants 



9 



x* and x'^. In order to do this matching we introduce the inner variable t in 
(5) defining it by 



t = 



(24) 



A 



Eq. (4) then rewrites 



u 2 - u tt - -rut = Q(x* + At). 



(25) 



Our nearest goal is to construct a part of the homoclinic (in terms of the slow 
variable x) orbit of Eq. (25) which tends to the solution (11) at large \t\ >> 1 
for | Xt |<S 1, and which is heteroclinic in fast variable t going from the saddle 
point to the stable focus (Fig. 2). Introducing the "energy integral" of Eq.(25) 



it is useful to evaluate Eq.(25) to the following integral for the function t(u). 
To the first order of approximation in <C 1 this function can be written 
down as 



Here Q* = Q(x*), = u(t = 0) and the integral must be taken along the 
orbit with changing the branches of the square root at the turning points, 
according to the sign of u t . The choice of the branches along the orbit can be 
understood with the help of Fig. 2. We have omitted the terms of order of A 
under the square root in (27) compared to /i/A by virtue of our assumption 
A 2 < ii. Iterating Eq.(27) we substitute u t from (26) at A, /i/A — > which 
yields for u t as zeroth order solution of Eq.(27): 



E — - 




= const, 



(26) 




(27) 




(28) 



Using normalized variable 




u 



(29) 
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Fig. 2. Phase plane of Eq.(25) 
instead of u and introducing 



from Eq.(27) we have as the second iteration (for negative t) 



<«) = -4/ dC 



Here 



2Q.V, \^(0 



where k* = k(x*) and C-i were defined by Eqs. (13,14) with 

3E 



2Q 



3 ■ 
2 



p(c) = c 3 -?c + ia+^g;5 S (c) 



In order to be matched with the solution (11) at t — > — oo, the solution ((t) of 
Eq.(31) must tend to the saddle point when t — > — oo (see Fig. 2). This means 
that both -P(C) and dP/d( must vanish at that point, which yields for the 
constant a 

a ~ 1 — — 

5QI A 



(33) 



Together with the expansions at k = k* ~ 1 

1 1 + kl 



C* = C-i(fc) = —4=77 - - 1 + l*?> ( 34 ) 



and 



• / I • / \ 27r\ 1 — a 

sm ( -arcsm(a) - — 1 ~ -1 + — — , (35) 



it defines the constant in Eq.(23) through the magnitude of k' 2 (x) at x = x*, 
i.e through k'^ = /^(x*), which is 



, /4 16V2^-i/i 



fc * =^-^ ^- (36) 



Under this choice of x'„ the asymptotics of the solutions (17) at £ — > — oo and 
(31) at £ — > — oo do coincide. Indeed, at large negative t the main contribution 
to the integral (31) comes from the vicinity of the point ( = 1/2. Expanding 
S(() at this point, so that 

^(o = (c-£) 2 (c-^ (37) 



where 

3 fi 



-1 + 



and substituting then a from (33) for large negative t we obtain 

V3 



Q*\J\~ 



: tanh 1 



\ 5 



(38) 
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Inverting this relation for u(t) yields 

(39) 

which is equivalent to the solution (11) or (17), taken at k ~ 1 and K£ in 
the argument of sn (see Eq.(17)) expanded at x ~ x* to the same order of 
approximation. 

Before we go to the extension of the solution to the neighbouring driver pe- 
riod, it is worthwhile to summarize the results of this section. Given the shock 
coordinate x* we obtained the solution (17, 23), which is valid throughout one 
period of the driver. Besides x*, the solution (17) depends on an additional 
constant x* (a point, where k(x) approaches unity and the adiabatic approxi- 
mation breaks down), which is connected with x* through Eqs.(22) and (36). 
This solution, however, does not allow the straightforward continuation to the 
nearby period of the driver, since at x = xo, where Qo = the adiabatic 
approximation breaks down again. This region needs to be treated separately 
which provides the matching conditions for the solutions of the form (17) in 
the two neighbouring driver periods. Therefore, this matching creates a link 
between the two subsequent shock coordinates x* and x*, which is the task of 
the next section. 



3 Matching the solutions 

In order to extend our solution beyond one driver period, let us turn now to 
the point x = xo, where the type of equilibrium between the nonlinearity and 
the driver changes from the stable to unstable one. Since u(x) tends to — \/Qo 
with increasing x — x*, we substitute 

u = R + v (40) 

into (5), where v(x) is supposed to be small, R = \^Qo denotes the regular 
branch of the square root at x = Xo which is positive at x > xo and is negative 
at x < xo (recall, that Qo(x) = Q ( x o) = 0). Then, linearizing (5) for small v, 
we obtain 

2Rv - X 2 v xx - fiv x = X 2 R XX + fiR x + Qi-\ , (41) 

where we denoted the constant 

Q 1 = Q - Qo ~ A 
13 



u ~ 2\/Q z + ( - — z } tanh z 




in Eq.(5). Introducing a new function 

= exp 



(42) 



we get the following equation 

4*xx 



2 i?+ ^ 2 



= F(x), 



(43) 



where Fix) denotes 



F(x) 



exp 



P 



-x 




2 Qi 



(44) 



Since both A and fx are small, one can show that to the leading order of 
approximation the solution of the Eq. (43) is defined by its homogeneous 
part, which can be written down asymptotically as 

<p(x) = a^j + S X , (45) 



where 



P 



1/4 



X = 



P 



1/4 



exp 



exp 



A 



S(x u x) 



A 



S(X!,X) 



(46) 



and 



S(x 1 ,x) = J v/prfx, 



p = 2R+^, p( Xl )=0. 



(47) 



Here X\ is the nearmost to x zero of the function p(x). The branches of 
y/p~ and constant c are chosen in such a way, that ImS* > | c |= 1 on 
each Stokes line k, originating at the point x\ and arg(cp -1 / 4 ) — > 0, x — > 
X\, x G k (Fig. 3). Note, that we defined the Stokes lines /j according to the 
condition ReS^i, x) = 0, x C Zj. Given the shock coordinate a;* the asymptotic 
solution is known to the left from x±, in particular on the line l\. A convenient 
representation of it can be derived from Eq.(17) for k 1: 
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Uq 



Fig. 3. The arrangement of the Stokes lines of Eq.(43) 



Here we have used the asymptotic solution of Eq.(22) for k 2 (x) in the region, 
where x satisfies the following two inequalities 



P(x - x'J > 1, x - x > exp —-(3(x - x[) 



(49) 



The latter restriction occurs since in the vicinity of the point x = x the 
solution k 2 (x) of Eq.(22) starts to grow becoming invalid because Qo(x) ap- 
proaches zero, as we noted earlier. 

To make the further analysis more transparent in addition to the inequalities 
(3,4) we introduce the restriction /x 2 A 3 , which is, however, not a critical one, 
and, at the same time corresponds to the most interesting case of developed 
oscillatory structures near the shock fronts. Comparing Eq.(48) with eqs.(40), 
(42), and (45), which represents the WKB asymptotics to the left from the 
turning point X\, we retrieve a and 8 in Eq.(45) 

a = -2 3 / 4 ^Ag^ /8 «) exp [Px'j2 - 2mP 1 /\] 

5 = a, (50) 
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where the bar denotes complex conjugate and 

X} Qt 

2 3 /2K(a;)(A; 2 + A;' 4 )V4- 



p /' dx f \ 

1 J 2V 2 K(x)(k 2 + k' 4 ) 1 / 4 ' 1 ; 



X* 



In our above calculations of the coefficients a and (3 we have extended the 
asymptotic solution (48) up to the point x — x\ (see Eq.(51)). If, however, 
xq — xi ^ exp —^(5{xq — x'J , which can be the case for relatively small 
x — x't (note that x — x\ ~ /i 2 /A 2 <C 1), the solution (48) together with the 
solution (11) becomes invalid and should be replaced by the solution which was 
obtained in this section for the vicinity of the turning point x — x\ (eqs.(40), 
(42), and (45)). According to this modification, the integrand in the phase 
integral (51) should be taken with k 2 (x) = in this region. This correction 
might become quantitatively important in calculation of the shock map for x* 
approaching xq (see ineq.(49)), which will be taken into account in the next 
section. 

The asymptotic solution to the right from x — x 1 can be obtained by making 
use of the standard WKB approximation, which links the coefficients a, 5 in 
the sectors D x and D 2 (see e.g. [16]). In the region D 2 the solution reads 

<f) = a 2 ^ + 5 2 x (52) 

with 

a 2 = e- iw/3 (a + ia) (53) 

and 

S 2 = ie~ iw/3 a. (54) 



Thus, a solution, which is valid also for real x > x\ can be written in terms 
of the function v(x) (Eq.(40)) as 

v{x) = 2- l ' i R- 1 l\x)e- px,2 ^ /i \(a + ia)e s{xi > x)/x + iae" 5 ^ 1 ^^] . (55) 



Taking fi 2 <^ A 3 into account, after a simple algebra from the last equation 
we obtain 



V57r Q!(x) 
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2cos( T P 1 + I jexp — 



JqI 



dx 



XI 



(in 2m „ \/2 f ± 

+ exp _ fl 



dx 



XI 



(56) 



Now we can construct the solution in the region x$ < x < x$ + 2n. Suppose 
the next shock is located at x = x* + 2tt, so that we take all the subsequent 
shock coordinates x (mod 2ir). Our task is to link x* with x*, which defines 
by an induction the solution for all x > x*. According to Sec. 2, the solution, 
which is valid in the region xo < x < xo + 2n, but not too close to the end 
points of this interval can be written in the form (17), with x* instead of x* in 
(21). Then, taking (11) at large negative £, i.e. using Eq.(39) with x* instead 
of x*, we obtain 



u(x) ~ y<5* < 1 — 12exp 



2 /l ; /sj X-X,- 

— x 



2tt 



(57) 



Here we introduced the notations = <3(x*) and ^ = C*(^*)> whereas the 
function is defined by Eq.(34). Expanding at small /i/A yields 



u ~ YQ* - 12 y exp 



1/4*^ 27T 



A 



P(x-x* -2n)/2 



• (58) 



This asymptotic solution should be matched with the asymptotics (56) in the 
region where they both are valid. This region is defined by the following two 
inequalities 

A x* + 2ir — x <C x* + 2n — x , 



To match Eqs.(56) and (58) we can evaluate the integral in Eq.(56) as follows 

X £*+27T 

J Ql /4 dx ~ J Ql /4 dx - Qi /4 (x* + 2vr - x). 

Xl Xl 



The asymptotic solutions (56) and (58) will coincide in the region of their 
validity if the following condition is met 



^ cos^P 1 + * 



V^Ql /8 VA 1 4 
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exp --(x* + 2vr - xj + — 



Ql^dx 



(59) 



XI 



4 Mapping the shock coordinate 



4-1 Arbitrary In- 'periodic driver 



Equation (59) provides an implicit relation between the two subsequent shock 
coordinates x* and x*. However, it is still rather complex for a comprehen- 
sive analysis. First of all there are two intermediate coordinates x' + and x\ in 
Eq.(59) that appeared during the evaluation of the map x* i— > x*. They have 
no other physical meaning than x* and x already do. Given the shock coordi- 
nate x*, the value x'^ can in principle be calculated by means of Eqs.(22) and 
(36), whereas x\ is defined by Eq.(47) and lies close to the point xo- Moreover, 
our assumption yU 2 <C A 3 allows us to replace x* by x* and x\ by xo, respec- 
tively (see Eqs.(22), (36) and (47)). We can also set xo = 2ir without loss of 
generality. After these simplifications we rewrite the map x* i— > x* as 

fM = #(x*), (60) 



where 



f(x) = Q 5 /8 (x) cos 



2tt 



TT 



Ql /4 dx 



+ 



TT 



,P(x/2-n) 



^2\J K{x){k 2 + A;' 4 ) 1 / 4 4 



(61) 



and 



9{x) 



\/30ttQ^ s (x) exp 



x/2 
A 



Q l J A dx + 



(3x 

T 



(62) 



4-2 Sinusoidal driver 



Let us consider the simplest case of sinusoidal driver in Eq.(l): 

Q' = \ sin(x) (63) 
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Fig. 4. a.) The functions f(x) and g{x) b.) The map of the shock coordinate F{x) 
Then, according to Eq. (6) and the text right after it we have for Q 

Qo = sin -, 

The functions / and g are shown in Fig. 4a. One can see, that given x*, Eq.( 60) 
can not be uniquely resolved for x*, since there are two different x for each 
g{x). At the same time, the function g(x) at x < x m (the point of maximum 
of function g(x)) corresponds to the position of the next shock x*, which 
is too close to xq to be correctly described by our approximation, since the 
asymptotic method of the derivation of the map x* \— > x* implies, that x* ^> A. 
Indeed, x m ~ A 2 / 3 and the solution with x* < x m can hardly satisfy the 
inequality x* A except for very small A. In general, this solution would 
correspond to the shockless one, which stays at the stable equilibrium u ~ 
— \[Q after the phase of the driver advanced beyond the point x (the jump 
would be of the order of x m — x <^ 1) and the matching should be treated 
on different basis. The above nonuniqueness is important when considering 
Eq.(5) per se. Here, we are looking for the travelling wave solutions of Eq.(l) 
and any particular choice of solution is justified. Interested in the shock-train 
solutions we choose G [x m ,27r] for <?(x*) from Eq.(60). In this case the 
mapping x* t— > x* is unique, although it is not defined for x, where 

/(x) < min (64) 

as well as in the interval [0,x m ). Resolving Eq.(60) numerically and omitting 
stars at x we get the mapping F(x) : x i— > x, which is shown in Fig. 4b. 
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5 Simple solutions of period one and the period-doubling bifurca- 
tions 



5.1 Weak chaos 



As discussed in the preceding section, the mapQ F(x) (Figs. 4) is not defined 
in those intervals where the inequality (64) is met. This means that after an 
iterate of the map F(x) has fallen into one of these intervals (denote them J^, 
and include also [0, x m ) in this category, as we have mentioned above), the 
travelling wave solution of the shock-train type cannot be continued to the 
right from this point. At the same time, a continuous to the zeroth order (in A 
and /i) solution could in principle be constructed in this driver period, as it was 
described in sect.l. However, we leave aside this possibility and restrict our 
consideration to the shock-trains. According to Eq.(60), the total number of 
intervals where F(x) is defined, can be estimatd as Ni ~ 1/A> 1. We denote 
these intervals as {Ik\k=i- The number of the fixed points of the map F{x) is 
N-p ^ 2N\. Each fixed point of F(x) corresponds to 2n -periodic solution of 
Eqs.(l) and (4). Since at a fixed point x 



ar J ax \= — — ~ exp 
g'(x) 



V2 
X 



x 

J Ql /A dx 



> 1, (65) 



these fixed points are mostly unstable with a possible exception for a point, 
where g(x) is not very small, or where f(x) reaches its maximum. Accord- 
ingly, the one parameter (x*) family of solutions of Eq.(5), constructed in the 
previous section with x* taken to be equal to one of this fixed points repre- 
sents an unstable periodic orbit of Eq.(5). It is important to note here, that 
the unstable in x periodic solution of Eq. (5) can correspond to the stable (in 
time) periodic travelling wave solution of the driven KdVB Eq.(l). This will 
be demonstrated in sec. 6. We also note, that the instability of a fixed point 
evolves in the positive x-direction, whereas the same fixed point attracts the 
backward iterations of the Poincare section even though the map has no unique 
inversion. The question how the spatial stability of our periodical solition is 
related to its temporal stability deserves a special treatment which is beyond 
the scope of this paper. Note, that even in linear systems the answer can be not 
a straightforward one (e.g. absolute versus convective instability in plasmas). 

We begin our analysis of the map F(x) with the case when a minimum of 
F(x) is close to one of its fixed points, say x = a, (Fig. 5). In this case the 
local behavior of the iterates near this fixed point does not differ from that 
which is known for one dimensional map with one maximum [17]. In particular, 

1 In what follows we relate x,F(x) to the interval [0, 1] instead of [0, 2tt]. 
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0.2 a 0.4 0.6 0.8 x 1 



Fig. 5. Mapping of the shock coordinate with one stable fixed point x = a for 
A = 0.27 and = 1.2. Horisontally aligned points, which are generated by F 100 (x), 
actually mark the basin of attraction of this fixed point 

starting with the stable fixed point x — a, | F'{a) |< 1, we obtain a series of 
period doubling bifurcations when, for instance, A is fixed A = 0.27 while j3 
changes from (3 = 1.2 down to (3 = 0.98 (see Fig.5-7). Below the point f3 = 0.98 
the local behavior of the mapping F(x) becomes chaotic. At the tangential 
bifurcation of the fixed point one can observe an intermittent behavior of the 
mapping [18], although the return time to the vicinity of the tangent point 
(where F(x) -i<Ci) might be quite large, since the orbit visits many other 
intervals of the mapping F(x) before it finds the return path, if ever. The 
measure one basin of attraction of the fixed point x = a and emerged 2 n 
cycles is made up of the points which are under the horizontal lines seen in 
Figs. 5-6 and under the scattered points shown in Figs. 7. When we start to 
iterate F(x), with an initial point from the basin of attraction the iterations 
will be attracted by the fixed point or will be trapped in 2 n cycle which is the 
fixed point of the map F n defined by F n = F o F n " 1 (x); F°(x) = x. Almost 
all other points of the [0, 1] interval will sooner or later escape in one of the 
intervals Jk where F(x) is not denned. This means that for almost all initial 
shock coordinates which are not from the basin of attraction, the shock-train 
solution can not be continued infinitely into the positive x-direction. The rest 
(a measure zero subset) of the [0, 1] interval, which we meant by saying "almost 
all" will be considered in the next section. It is worth noting, that although 
the measure zero solutions are exceptional for the ODE (5), they can be very 
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0.2 a 0.4 0.6 0.8 1 

X 



Fig. 6. Stable cycle of period 2 near the fixed point a, X = 0.27, (3 = 1.1. The 100-th 
(F wo (x)) iteration is shown as well as F(x) 



a.) b.) 




Fig. 7. a.) The same as in Fig.6 but for (3 = 0.98 b.) A magnification of the box in 
Fig. 7a. 

persistent and quite "typical" travelling wave solutions of the correspondent 
evolution equation Eq.(l). 
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Fig. 8. A piece of the map F(x) illustrating the fragmentation of phase space under 
iterations. The points from the black areas will escape as it is shown by the arrows. 



5.2 Strong spatial chaos 



The behavior of the iterates considered in the previous section is chaotic for 
j3 below f3 c = 0.98, but this chaos is weak in the sense that the spatial profile 
of the solution of Eq.(l), corresponding to this sequence of shock coordinates 
{x n }^ =0 is still rather close to periodic solution since | x n — a A. By further 
decreasing of (3, the measure one basin of this weakly chaotic attractor around 
x = a disappears, and almost all initial points will escape in [0, 1] \ \Jk=i Ik 
as it was argued in the previous subsection. To identify the nonwandering 
points of F which define now the sequence of shocks in the spatially chaotic 
solution of Eq.(l), we proceed as follows. Consider the nonescape set, which 
occurs under the backward iterations of U^i h an d denote its limit by A = 
[0,1] \ U^Lo F~ n (1J (recall, that J k denotes the escape intervals). The 
construction of this set is similar to that of the classical Cantor set and is 
illustrated in Fig. 8. Once the initial shock coordinate is taken from the set A, 
the shock train solution can be continued ad infinitum. 

An important property of this set is its capacity dimension. One calculates it 
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Fig. 9. Illustration of the natural measure of chaotic repellor of the map F(x). All 
four plots show the lengths of intervals that have not escaped after six iterates of F, 
as a function of their left boundary. Each box is magnified in the subsequent plot. 



from the following formula, e.g. [19,20] 



lim logiV£ 
~o log(l/e) 



(66) 



where A" is the minimal number of the intervals of length e, needed to cover 
the set. To calculate D c we approximate the set A by iterating the interval 
[0, 1] and by removing all those x, which after n iterations leave the domain 
of F, i.e. escape in one of the intervals Jk- Denote the rest of [0, 1] by A n = 
[0, 1] \ Urn= F- m (U Jk) = {(xi,Xi + A Xl )}^ v so that A = lim^oe A n . The 
distribution of segment lengths Axi (N s is their total number) of the set 
remained after six iterations Aq is shown in Fig. 9 as a function of the left 
boundary of each segment. The iterations A n cannot be continued beyond 
certain n = n max , because of the limited computation accuracy. Since the map 
F(x) expands phase space very quickly (see Eq.(65)), errors in the calculation 
of the points from invariant set of F will rapidly grow and after a few iterations 
the point escapes into one of the intervals Jk just due to insufficient accuracy. 
Therefore, for numerical purposes we consider the sets A n with n = 5, 6, 7. 
After n ^ 8 iterates the number of segments N s in A n starts to decrease for 
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Fig. 10. a.) The number of intervals needed to cover the sets A n for n = 5, 6, 7 and 
for A = 0.3, = 1.3 b.) A linear fit to lniV for the case n = 6. 

above reasons. The dependence N(e) for different A n is shown in Fig. 10. Since 
A n is measure one for any finite n , even though very irregular, the dependence 
iV(e) becomes N ~ e~ l for e smaller than typical segment length in A n . At 
larger e one can see a distinct region of the linear behavior in the log — log 
scale, which captures the main fractal properties of the limiting (n — > oo) 
set A. In this e- range the slopes of the curves are essentially identical for 
4 < n < 8. The transition to the D c = 1 behavior at smaller e occurs earlier 
for smaller n which is a natural result. For n > 8 the slope becomes smaller, 
which underestimats the capacity dimension of set A due to the escape of 
points from A n because of the limited computation accuracy. For A = 0.3 and 
fj — 1.3 the capacity dimension D c m 0.26. 



6 Comparison with the time dependent numerical solutions 



Before we verify numerically some of our above results it is worthwhile to 
comment on the possible approaches to this verification. First of all there are 
two different aspects of our analytic calculations. One of them concerns the 
response of damped nonlinear oscillator (5) to the large amplitude forcing. 
In the present paper we leave aside the problem of numerical study of (5) 
regarding its reduction to the discrete map (60). We only note here that such 
a study would not be straightforward, since the typical fixed point of this 
map is unstable as well as the motion on the repellor, considered in the pre- 
vious section and the solutions of (5) are unbounded in their vast majority. 
Backward iterations of this map, locally attracted by this fixed points, would 
require a choice of preimages, i.e. the branches of F~ 1 (x), since F(x) is not 
invertible. Even though this numerical study were successively done, it would 
not give us an ultimate answer to an important question concerning the driven 
KdVB equation. In particular, it would be still unclear which solutions of the 
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travelling wave type, i.e. the solutions of (5) are the attractors of Eq.(l)? We 
therefore perform direct integration of the evolution equation (1) in time. It 
should be realized, however, that this approach is also insufficient for a firm 
manifestation of spatial chaos, since it would require an infinite number of 
modes or equivalently an infinitely long system. An additional problem is that 
we already have two essentially different scales (dispersion/viscosity and the 
driver scales) and in order to detect at least a possible trend to the spatially 
chaotic behavior we would need to model a system that appreciably exceeds 
the driver scale. Note, that it is the shock-train technique, developed in the 
previous sections that is especially suitable for handling systems with consider- 
ably different spatial scales. At the same time, the spatially chaotic travelling 
wave solutions, obtained in the preceeding sections are not necesseraly the at- 
tractors of time dependent Eq.(l). An example of the impossibility of chaotic 
solutions is given in [21] where the Frenkel-Kontorova model is considered. In 
particular, it was proven in [21], that the ground-states of this model can be 
periodic or quasi-periodic but the spatial chaos is ruled out. 

The numerical investigation of the onset of spatial chaos is beyond the scope 
of this paper. The above arguments suggest that it could be a rather difficult 
task and certainly deserves a separate study. First of all one should check 
whether at least the first steps in a route to spatial chaos can be described 
correctly by an analytical shock-train solution obtained in previous sections. A 
possible strategy consists of obtaining 2n periodic solutions numericallly and 
of further identification them with the correspondent fixed points of the map 
F(x) (sec. 4). Then, the 4-7r-periodic solutions should be found and identified 
with the fixed points of F 2 = FoF. These fixed points appear in the vicinities 
of the unstable fixed points of F. A continuation of this process should in 
principle give rise to the spatially chaotic solutions. This section deals with 
the very begining of this program. First, we demonstrate that the travelling 
wave solutions of (5) do attract time dependent solutions of Eq.(l). Second, 
an approximate Poincare map can really capture their complexity. 

We start with the numerical integration of Eq.(l) for periodic driver given 
by Eq.(63) with different initial conditions. The results are shown in Figs. 11. 
One can see from these two figures that the shock-train is a persistent solu- 
tion of Eq.(l). Starting from large amplitude sinusoidal profile (Fig. 11a), the 
solution steepens very quickly and forms a shock-train which has the driver 
period. Fig. lib demonstrates the case with completely different initial con- 
dition, where a rather random small amplitude perturbation is imposed. It 
decays rapidly due to the dispersive spreading and viscosity. Then after the 
amplitude of the driver's harmonic achieved nonlinear level, it steepens to 
form exactly the same travelling wave solution. Figs. 12 show the compari- 
son of such a steady state with an asymptotic solution given by the formulae 
(17,21-23). It is seen from Fig.l2a that the agreement is very good, taking 
into account that only zero order solution u (£,x) is plotted. There are, how- 
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Fig. 11. Formation of the shock-train for two different initial conditions, and the 
driver given by Eq.(63) The coordinate is shifted back by ir/2 and compressed to 
[0, 1]. Only one period of the solution is shown. 



ever, two regions where the deviation is significant. They correspond to the 
nonuniformities of the asymptotic solution (17), resolved in sees. 2 and 3, but 
not taken into account in Figs. 12 for simplicity. First, the amplitude of the 
first peak is somewhat larger than numerical one due to our replacement x'^ 
by x*, which overestimates k(x) at x ~ x* (see sec. 2). Second, the deviation 
in the region where an oscillatory behavior is changed for a monotonic one 
results simply from invalidity of Uo(£,x) given by formula (17) at the turning 
point x (sect. 3) in the case of moderate (3. (We have taken (3 = 1.8 to obtain 
pronounced solitons between the shock and the turning point x ). Thus, the 
agreement can be even improved by using the asymptotics given in sec. 2, 3 
instead of Eq.(17). 

The solution uq(£,,x) depends also on additional parameter, the shock coor- 
dinate x*. For the period-one solutions considered here, x* must be equal to 
one of the fixed points of the map F(x), plotted in Fig. 12b. In the analytical 
solution shown in Fig. 12a, x* = 0.497 which corresponds to the filled circle in 
Fig. 12b. Solutions corresponding to the empty circles in Fig. 12b are shown in 
Fig. 12c and they clearly do not agree with the numerical solution. This allows 
us to identify a continuous solution with one of the discrete shock coordinates. 
There are as many as twenty fixed points of the map F(x) shown in Fig. 12b. 
Each of them corresponds to a period-one solution of the driven KdVB equa- 
tion. Moreover, one can easily draw F n (x) for n > 1 and obtain, in principle, 
all the solutions of period n, as the fixed points of F n . Due to the limitations 
of applied asymptotic method, a few right most fixed points cannot be de- 
scribed correctly. At the same time they correspond to the solution, in which 
u ~ yfQ throughout the driver period, i.e. actually to the counterpart of the 
solution u ~ —\[Q briefly discussed in sec. 4. 2. This situation can in principle 
be systematically improved by adopting an appropriate asymptotic method 
for these solutions, which would allow us to incorporate them in to the shock 
train class simply as limiting solutions with the vanishing shock amplitude. 

The fixed points of the map F(x) produce a discrete set of 27r-periodic so- 
lutions. Each solution contains an integral number of solitons that gradually 
decrease in amplitude filling up the region between the shock front and the 
turning point. Why does the system follow the solution with x* = 0.497? 
It is clear that the system choose a steady state to follow according to the 
constraints imposed by the initial conditions. Such a constraint is the mass 
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a.) 




Fig. 12. a.) Zeroth order analytic solution Eq.(17) (solid line) having the shock at 
= 0.497 (the correspondent fixed point of the map F(x) is marked by the filled 
circle in Fig. 12b). Dashed line shows the steady state numerical solution. Both 
curves are drawn for the case A = 0.1 and = 1.8. 

b.) The map F(x) for the same A and (3. c.) The same numerical solution as in 
Fig. 12a (heavy line) and the two analytical solutions (thin and dotted lines) having 
shocks at x* marked by the opened circles in Fig. 12b. 

invariant 

i 

u = Ju(x,t = 0)dx, (67) 
o 

which was set to zero. In order to estimate u from the lowest order solution 
(17), we can to the good approximation average the fast oscillations and regard 
k as a small value for x > x*. The solution u (£(x*,x),x) will then produce 
u = only for x* = 1/2. For arbitrary u- value, the system follows the solution 
with an appropriate x*-value. It is interesting to note that a solution taken 
from the set u (£(x n , x), x), being discrete fulfills the condition u = only 
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approximately, although with a good accuracy ~ A 1 (here x n is one of the 
fixed points of F(x)). Thus, the difference might in principle be compensated 
by the higher order terms. We do not consider this problem here. More detailed 
comparison of our analytical solutions with the numerical ones including the 
doublings of the spatial period will be done in a future paper. 



7 Discussion 

An interesting feature of the repellor considered in sec. 5. 2 consists in the 
concentration of its natural measure inside narrow intervals where F'(x) ^> 1. 
The function F(x) maps these intervals onto the appreciable part of [0, 1] and, 
therefore, those are the points from these intervals which have better chances 
to avoid the escape set J k after a number of iterates. All the unstable fixed 
points, shown in Fig. 12b belong to these intervals. Motion on the repellor is 
essentially a random walk between these narrow intervals. If we neglect their 
width together with the underlying multifractal structure of these intervals, we 
obtain a discretization of the shock coordinates in the shock-train. The mean- 
ing of such a discretization is obvious; it corresponds to an integral number 
of vibrations or solitons between the shock coordinate and the turning point 
x . These regions of oscillatory behaviour can also be considered as optically 
transparent domains, created nonlinearly by the driver ("potential wells"). 
This can be easily understood from the consideration of the solution near the 
turning point x , where the small amplitude oscillations allow the linear treat- 
ment (sec. 3). Figs. 12b, c in particular demonstrate that at each second fixed 
point of F(x) the number of solitons differs by unity. Note, that the majority 
of the fixed points are close to zeroes of function f(x) in Eq.(60). This fact 
suggests a simplified "quantization rule" for the discrete shock coordinate x* 
in the shock-train (which is valid for sufficiently small g(x)) 



In general, even when the parameter (3 is not large and the linear approxi- 
mation at the turning point xo is not justified, the number of solitons is a 
constitutive aspect of the considered chaos scenario. To put it anotherway, 
consider the motion on the invariant set of the map F(x). Suppose we know 
the coordinate x* of some reference shock in the shock-train or, equivalently, 
the coordinate of the leftmost soliton in the correspondent potential well. An 
appropriate accuracy to which x* should be known is o(A), as Eq.(68) sug- 
gests. In practice, that would be enough only for identifying x* with one of 
the region of the concentration of the natural measure, marked by the digits 
in Fig. 9a. Better accuracy in shock coordinate would allow us to go only few 





29 



steps deeper into coarsegraining of this region. Since the Lyapunov exponents 
are typically of order of 1/ A on the repellor, we cannot even predict into which 
region falls the coordinate of the next shock in the shock-train. With a more 
accurate knowledge of the reference shock coordinate we could at best pre- 
dict a few next. Therefore we will obtain an apparently random sequence of 
soliton numbers in the subsequent potential wells, created by fully determin- 
istic driver. In general, one can argue in terms of the topological conjugacy 
between these random sequences and the set of real numbers, as it is widely 
used in study of deterministic chaos in quadratic or binary tent maps, see e.g. 
[22,23]. The number of the coarsegrained intervals in Fig. 9a or, equivalently, 
the number of solutions of Eq.(68) for be taken as a base of arithmetic. 



8 Conclusions 

Using asymptotic methods we have represented the continuous travelling wave 
solutions of the periodically driven KdV-Burgers equation as a Poincare map 
at the driver period. Unlike many other analytical approaches to driven equa- 
tions, the amplitude of the driver is not assumed to be small. We argue that 
the typical response of this system to the external force consists of a frag- 
mentation of the solution into two alternating parts. One of them is smooth 
and has the driver scale. The second is an oscillatory one (shock), its scale 
is defined by the small dispersion and viscosity. The scenario of the spatial 
chaotization of this solution is suggested. It consists in a randomization of 
the shock-coordinate with respect to the driver phase under the action of the 
derived Poincare map. 
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